function [wdraw,what_star,wplus]=kalman_dr1(y,xdata,Z,Cobs,T,R,Qmat,A,addin)
% ==================================================================
% KALMAN_DR1
% 
% This function takes as inputs the matrices of a state space model with no
% error in the observation equation and the Kalman Gain, and F^(-1) matrices from 
% a forward run of the Kalman Filter in order to simulate a 1 single set of draws
% from the state disturbances and the states using a disturbance smoother. 
% It also returns the filetered and smoothed states. 
% It works both with ESTDGSE and SSM type models. 
% Only for NON_DIFFUSE state space with constant variances 
%
% The models is 
% y(t) = Z*s(t) 
% s(t) = T*s(t) + A*xdata(t) + R*n(t)  V( n(t) ) = Q=Qmat'*Qmat   
%
% Where y is [NZ  1], xdata(t) [NE  1]; s(t)[NS 1] and n(t) [NX x 1] 
% 
% Related files: 
% KALMAN_DRALL.m  Generate multiple draws 
% Inputs 
% ======
% y     nobs x NZ stacked vector of observations, nobs is the total number of observations  
% 
% All matrices but Q(t) are assumed to be time invariant to simplify the computation 
% The correct size of input Qmat then depends on whether Q is time varyong or  not. 
% not time varying    Qmat is NX x NX 
%      time varying   Qmat is NX x NX x T where
% 
% Alejandro Justiniano and Giorgio Primiceri (c) April 2005 
% =============================================================

% 1. Initial check for dimension
% -------------------------------
[nobs,nz]=size(y);
[ny,nx] = size(R);
ns=size(T,1);
Q=Qmat'*Qmat;

% 2. Removal of Auxiliary data and means
% ydem
% ---------------------------------------
if ~isempty(xdata);
    flag_xdata=1;
    check   =size(xdata);
    if check(1)~=nobs;error('X and Y must have same number of rows');end
    if check(2)~=size(A,2);error('Cols A must match Cols of X');end;
    if size(A,1)~=size(T,1);error('Rows A must match dimension state space');end;
    ydem=y-xdata*(A');
else
    flag_xdata=0;
    ydem=y;
end

if any(Cobs~=0)==1;
    ydem=ydem-repmat((Z*Cobs)',[nobs 1]);
end;

% 3. Initialization
% ==================
% I. Non-diffuse
% PSHAT and SHAT
if addin.flag_diffuse==0
    if addin.flag_pzero==1
        pzero=addin.pzero;
    else
        pzero=disclyap_fast(T,R*Q*(R'));
        if ~isempty(find(isnan(pzero))~=0)
            error('PSHAT contains Nan entries')
        end
    end
else
    disp('Still have not added diffuse filter')
end
if addin.flag_szero==1;
    szero=addin.szero;
    if length(shat)~=ns;error('SHAT must be NSx1');end
else
    szero=zeros(ns,1);
end

% ==============================================
% 4.    Simulate the innovations and the data
% w+
% y+
% s+    alpha+
% =============================================
Ztr=Z'; %clear Z;
Rtr=R'; %clear R;
Ttr=T'; %clear T;


wplus=randn(nobs,nx)*Qmat;
% This is y+ and alpha+ in DK.
% Note that the first state is drawn from a Normal (shat,pshat)
[yplus,stplus]=simssmodel(wplus,T,R,Z,szero,pzero);

clear rstar;
% Forward with ydem
[what,shat]=forward_backward(ydem,szero,pzero);


clear rstar;
% Forward with yplus
[whatplus,shatplus]=forward_backward(yplus,szero,pzero);
shatplus=shatplus';

plot([shatplus(:,1) stplus(:,1)])
% According to the paper, if a single draw is needed,
% ystar=y-yplus generated above
wdr=what-whatplus+wplus;
sdr=shat-shatplus+splus;

    function [wsmooth,asmooth]=forward_backward(yy,a_zero,p_zero)
        wsmooth=zeros(nx,nobs);
        asmooth=zeros(ns,nobs);
        % Forward recursion
        %[junk,vstar,finmat,kmat]=kfilter_mod(yy,Z,T,R,sqrt(2)*Qmat);
        [junk,vstar,finmat,kmat]=kfilter_mod(yy,Z,T,R,Qmat);
        % Backward
        
        rstar=zeros(ns,1);
        ii=nobs;
        for ii=nobs:-1:1;
            Ltr=Ttr-Ztr*kmat(:,:,ii)';
            [rstar,wsmooth(:,ii)] = smoothdis( rstar, Q, Rtr, Ztr, finmat(:,:,ii), Ltr , vstar(:,ii) );
        end
        % Forward state
        asmooth(:,1)=T*a_zero+p_zero*rstar;
        ycheck(:,1 )=Z*asmooth(:,1);
        
        ycheck=zeros(nz,nobs);
        for ii=2:nobs;
            asmooth(:,ii)=T*asmooth(:,ii-1)+R*(wsmooth(:,ii));
            ycheck(:,ii) =Z*asmooth(:,ii);
        end
        
        check=max( abs(ycheck'-yy)  )
    end
end